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ABSTRACT 

We examine the nonlinear development of unstable magnetosonic waves driven by a background radiative 
flux - the Radiation-Driven Magneto-Acoustic Instability (RMI, a.k.a. the "photon bubble" instability). The 
RMI may serve as a persistent source of density, radiative flux, and magnetic field fluctuations in stably- 
stratified, optically-thick media. The conditions for instability are present in a variety of astrophysical environ- 
ments, and do not require the radiation pressure to dominate or the magnetic field to be strong. Here we nu- 
merically study the saturation properties of the RMI, covering three orders of magnitude in the relative strength 
of radiation, magnetic field, and gas energies. Two-dimensional, time-dependent radiation-MHD simulations 
of local, stably-stratified domains are conducted with Zeus-MP in the optically-thick, highly-conducting limit. 
Our results confirm the theoretical expectations of Blaes and Socrates (2003) in that the RMI operates even 
in gas pressure-dominated environments that are weakly magnetized. The saturation amplitude is a monotoni- 
cally increasing function of the ratio of radiation to gas pressure. Keeping this ratio constant, we find that the 
saturation amplitude peaks when the magnetic pressure is comparable to the radiation pressure. We discuss the 
implications of our results for the dynamics of magnetized stellar envelopes, where the RMI should act as a 
source of sub-photospheric perturbations. 

Subject headings: radiative transfer — diffusion — MHD — instabilities — methods: numerical 



1. INTRODUCTION 

Many gravitationally-bound astrophysical systems remove 
their binding energy by the diffusion of radiation through op- 
tically thick regions. Examples are main-sequence stars and 
accretion disks onto compact objects. In general, these sys- 
tems tend to be good conductors and therefore support mag- 
n etic fields. 

iBlaes & Socratesl (12003b hereafter BS03) found that these 
environments are susceptible to magnetosonic overstability. 
The physical mechanism responsible for secular driving in- 
volves short wavelength compressible fluctuations that grow 
exponentially due to the presence of a background radiative 
flux. Fluid motion along the direction of the equilibrium mag- 
netic field, resulting from magnetic tension forces, couples to 
the radiative flux perturbation. In the event that these two per- 
turbations are in phase, the radiation field performs work on 
the fluid oscillation and in crease s its amplitude (Figure[T}. 

A surprising result of BS03 is that weakly magnetized 
and/or gas pressure dominated equilibria are subject to the 
same radiative driving mechanism responsible for oversta- 
bility in radiation pressure and/or magnetic pressure domi- 
nated environments. The instability mechanism is so generic 
that even when the radiation is degener ate (e.g., neutrinos i n 
core-collapse supernova environments. ISocrates et aI1l2005l) . 
it may still operate. Furthermore, a similar mechanism 
operates in the effect that the energy-transporting particles 
are charged and diff use primarily along magnetic field lines 
(ISocrates et al .1 1 20081) . 

Before lBS03l such phenomena were thought to be restricted 
to radiation-pressure-dominated media that are strongly mag- 
netized. For these conditi ons, the instability is often referred 
to as " photon bu bbles" (iPrendergast & Spie sell 1 1 97 3t lArons I 
119921: iGammid [1998; Blaes & Socratesl 1200 U iBegelmanl 
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120011) . though neither buoyancy, nucleation, nor surface ten- 
sion play a role in the driving mechanism. Numerical in- 
vestigation of the nonlinear development of this instability 
has been limited, with only a few calculations of strongly- 
magnetized and mostly radiatio n pressure dominated atmo 
spheres reported in the literatu re ( Klein gtinil996HHsu et 
1997; Davis et al.|2004t|Turner et alj2004 120051: iTurner et 
2007; [Jiang et alj 120121) . A major difficulty for numerical 
studies results from the fact that the radiation diffusion time 
over a wavelength smaller than the gas pressure scale height 
needs to be resolve d in order to attain numerical convergence 
(ITurner et al.l 120051) . In practice, the oscillation period itself 
is often much shorter than the e-folding time of the mode, re- 
quiring a large number of time steps to capture the secular 
amplification of the waves. 

Here we take advantage of the parallel radiation-MHD code 
Zeus-MP (Hav es et al.l 120 06) to begin an exploration of the 
more extended parameter space for unstable modes found by 
IBS03L To reflect the physical character of the driving mecha- 
nism, we designate it the Radiation-Driven Magneto-Acoustic 
Instability (RMI). 

The parameter space to be explored is motivated by the ra- 
diative envelopes of massive stars, many of whic h are known 
to sup port up to kG-strength magnetic fields (e.g., Wa de et all 
2012 and references therein). In main-sequence O stars, the 
radiation pressure only marginally exceeds the equilibrium 
gas pressure; a magnetic field near equipartition with the gas 
provides favorable condit ions for the growth of the RMI (e.g., 
IBS03L ITurner et ai]|2004l) . Going down the main sequence in 
mass and luminosity, radiation pressure becomes increasingly 
less important, yet the instability may persist near the photo- 
sphere. 

The paper is organized as follows. Section 2 contains an 
overview of the physics of the RMI. Section 3 describes our 
assumptions, the numerical method employed, and the param- 
eter space explored. Section 4 presents our results. Section 5 
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FIG. 1. — Driving of magnetosonic modes by a background radiative flux 
F (BS03). The shaded region represents a magnetosonic plane wave, with 
darker and lighter regions corresponding to overdensities and rarefactions, 
respectively. The perturbation to the flux <5F does work on the component of 
the velocity perturbation <5v perpendicular to the wave vector k, provided that 
the latter is not entirely parallel or perpendicular to the magnetic field B. The 
acceleration of gravity is denoted by g, and velocity vectors are drawn for a 
slow mode in the weak field limit. 



briefly discusses implications for massive stars. We conclude 
with a summary of our work in Section 6. Appendix [A] con- 
tains a brief derivation of the growth rate of the RMI and the 
stability criteria. Appendix IE1 describes verification tests of 
Zeus-MP. Mathematical symbols have their usual meanings; 
the most frequently used ones are defined in Table Q] 

2. OVERVIEW OF THE RMI 

The photon bu bbl e instabili t y wa s discovered numerically 
by lAronsl (fl992) and lGammid (119981) in atmospheres that are 
radiation pressure dominated, with either strong (B 2 /8tt 3> 
p g ) or super-strong (B 2 /tt ^ p g + E/3) magnetic fields. In 
their computations, they found the slow magnetosonic wave 
to be susceptible to over-stability, driven by the p r esence 
of a background radiative flux F. iBlaes & Socratesl (120011) 
subsequently found that the fast magnetosonic mode could 
also be driven unstable, and that in the presence of rota- 
tion, the splitting between the Alfven and slow wave al- 
lows for the - otherwise incompressible - Alfven branch to 
be driven as well. The effects of magnetic stratification in 
the magnet ic- and radia t ion-do minated regime have been ad- 
dressed by lTao & B laes (20111). 

The stability criteria and dynamics of the RMI for a con- 
stant magnetic field were derived in BS03, They found that (i) 
stability is determined by a balance between driving from the 
background flux F, damping due to radiative diffusion, and 
departures from thermal equilibrium, (ii) the instability can 
be thought of as magnetosonic waves that are secularly driven 
by the background flux F, (iii) the instability mechanism oper- 
ates independently of the degree of thermal coupling between 
the gas and radiation, and (iv) the instability mechanism oper- 
ates even when the magnetic and/or radiation pressure is/are 
less than the equilibrium gas pressure. 

In order to determine the relevant parameter space for nu- 
merical simulations of the RMI, we must understand the 
physics that determines whether a given atmosphere is unsta- 
ble to RMI driving. A brief derivation is found in AppendixlAl 
Below we provide a general overview. 



TABLE 1 
Frequently-Used Symbols 



Symbol 


Quantity 


P 


density 


v 


velocity 


e 


gas energy density 


E 


radiation energy density 


T 


temperature 


s 


entropy per baryon 


n 


baryon number density 


Pg 


gas pressure 


Ci 


isothermal gas sound speed 


g 


gravitational acceleration 


Pt 


radiation pressure 


F 


radiative flux 




flux-mean opacity 




absorption opacity 


B 


magnetic field 


V A 


Alfven velocity 


Pm 


magnetic pressure 


b 
ft 


wave vector 


U! 


mode frequency 


^diff 


diffusion frequency (eq. (2)) 


^Ih 


heat exchange frequency (eq. (4)) 




w 2 -(k-v A ) 2 


3i 


rapid-diffusion parameter (eq. 1141) 


Sc 


critical Eddington ratio (eq. 1151) 


c 


speed of light 


L 


box size 


CsO 


adiabatic gas sound speed at z = 


Po 


density at z = 



2.1. Characteristic Frequencies 

We consider here the evolution of plane waves of frequency 
u) and wave vector k on an ideal MHD fluid with constant 
background magnetic field B and gravitational acceleration g. 
Furthermore, we take the various frequency-weighted opaci- 
ties to be constant in order to isolate the RMI from opacity- 
driven instabilities. 

All of the unstable RMI parameter regimes found by BS03 
involve heat rapidly diffusing across a wavelength in compar- 
ison to the frequency of oscillation. That is 

Wdiff > Re(w), (1) 

where 

^diff = 5 (2) 

is the radiative diffusion frequency, and Kf is the flux-mean 
opacity. Radiation diffuses rapidly in comparison to the 
mode oscillation time, i.e., the perturbations are highly non- 
adiabatic. 

In this paper, we restrict ourselves to the case where the 
gas and radiation are thermally well coupled, so that both are 
described by a single temperature. This condition is satisfied 
in the radiative envelopes of massive stars. Quantitatively, we 
demand that 

Wb, » Re(w), (3) 

where 

4(7- l)g 

W t h = KaPC, (4) 

Pg 

is a characteristic heat exchange frequency (BS03), with n a 
the thermal absorption (Planck-mean) opacity, and c the speed 
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of light. In this limit, gas-acoustic perturbations travel at the 
isothermal sound speed, cf = p g /p. 

2.2. Stability Criteria 

As shown in Appendix lAl the criterion for RMI driving is 
given by (eqns. IA201 - 1A2T1 ) 

[ W 2 -(k-v A ) 2 ] 



kF > LonT 



ds 



dlnp 



(5) 



The expression above quantifies the competition between the 
rate of radiative energy transfer into a fluctuation over its 
wavelength \ = 2tt fk due to the presence of a background flux 
F, and the rapid diffusion of energy that results from com- 
pression. Equation (O is obtained as a first order expansion in 
the small parameters (kH g ) , Im(w)/Re(w), and Re(w)/o; t iiff, 
where H g is the gas p ressure scale height. Equation (O differs 
from equation ( IA201 > only by angular factors, being thus ac- 
curate to withi n a fac tor of order unity. A more accurate value 
than equation ( IA201 > can b e obta ined by solving the eighth- 
order dispersion relation of lBS03l 
For slow magnetosonic waves, we may write 



4 

p+-E 



Ci, 



where we have used equation (IA16b and defined 



( = min I — , 1 

Ci 



The growth rate can be approximated by 



where 



4c ; V 4E J s 

kfF 

e = 

eg 



F 



F 

Fedd 



p+\e 



(6) 



(7) 



(8) 



(9) 



is th e Edd ington ratio. The angular dependence of equa- 
tion dA20b has been replaced by a global factor of 1/2 that 
accounts for the product of the sines of the relative angles be- 
tween k, B, and F. 

For fast magnetosonic waves, the corresponding instability 
criterion reads 



F> 



1 



4 

P+~E 



(10) 



with a corresponding asymptotic growth rate 



Tfast : 



EL 



1 + 



3p 
4E 



£1 

F 



4 

P + - 3 E 



va 



(11) 

The extra factor of 1/2 has also been included. 

It is instructive to visualize the instability criteria for slow 
and fast waves in a plane with the Eddington ratio e on the 
x-a xis a nd the ratio va/c, on the y-axis, as is shown in Fig- 
ure [2/2] The critical stability curves can be obtained by relat- 
ing the Eddington ratio to the ratio of radiation to gas pressure 
through 

4,52-^. (.2) 
p l-e 

3 This is not the same as using an isothermal equation of state, since oj^ is 
a local quantity. 




£ = F/F 



FIG. 2. — Unstable parameter regions in the thermally locked case (mode 
frequ enci es smaller than eq. (4])- If the "rapid-diffusion" parameter 5ft 
(eq. 1141 ) exceeds unity, then a critical value of the Eddington ratio e c 
(eq. 1151 ) exists. For e > e c , unstable fast-magnetosonic modes are possible, 
and slow-magnetosonic modes are unstable for any magnetic field. Below e c , 
unstable slow modes are possible only up to a maximum value of v A / Cj given 
by equation (6), and fast modes are stable. If 5ft < 1, then for any Eddington 
ratio there exists a maximum value of v A fa < 1 for which slow modes can 
be destabilized. 



where 



3l 

H r 



(13) 



is the ratio of gas pressure scale height to radiation pressure 
scale height H r . 
The resulting curves have an additional free parameter, 



k>- 1 JL c JE± 
4 r ' 



(14) 



where t p = KpH g is the optical depth over a gas pressure scale 
height. The parameter is proportional to the ratio of diffu- 
sion speed c/t p to the isothermal sound speed. When > 1, 
there is a critical value of the Eddington ratio 



£c = 



Tp Ih 



(15) 



above which fast and slow waves are unstable for any mag- 
netic field strength. If 5ft < 1, then only slow modes can be 
destabilized, and for every Eddington ratio there is a maxi- 
mum field strength for which the slow wave RMI operates. 

2.3. Aspects of the Stability Criteria 

The growth/damping rate and stability criteria above give 
insight regarding the properties of the equilibrium about 
which the RMI takes place. Balance between radiative driv- 
ing and damping is, respectively, a balance between energy 
input from the background radiative flux F and energy loss 
due to diffusive cooling that results from compression, which 
is oc (ds/dp)r- 

Consider, for example, the slow magnetosonic wave in the 
strong field limit, so that equation <(5j becomes 



1 ^ 1 



1 Ec^ 1 , , 



ds \ 
dlnp J T 



4 
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> 



At 



diff 



(Pg + 4p r ) 
At u , 



(16) 



dE 4 

— + v-VE= — EV-\-n R pc (E-aT?)-V -F (21) 
ot 3 s 



where At = 27r/A; is the wavelength of an oscillation with wave 
number fc, and At u = 2tt/lu is the oscillation period. The 
left hand side is the upper limit - if, e.g., the perturbations 
were nonlinear - for the rate of radiative energy transfer into 
a fluctuation of scale A^.. The characteristic velocity at which 
the background radiative field transfers radiative energy from 
one layer to another is given by c/r p , and the corresponding 
time at which radiative diffusion transfers energy across the 
wavelength A^ is given by Afdiff. The right hand side of the 
criteria for instability results from entropy or pressure pertur- 
bations that are sourced by changes in density. Since radia- 
tion diffuses, effectively, instantaneously over a wavelength 
Aj during the oscillation time At u , the rate of energy loss 
per unit volume is capped by the relative contribution to pres- 
sure or entropy changes that arise from compression, which is 
oc (ds/dp)r and the oscillation time At u . 

By re-inserting S (eq. IIA22I ). the instability criteria can be 
written as 

-r— >-s[w,k,B] 

Af diff £ 



^>-H[,k,B] 



(17) 



where the phase velocity v p h = u/k (eq. IA4II ). and we have 
taken ep g /p r ~ 1. The maximum value of the left hand side 
is achieved at the photosphere where r = 1 . 

Equation (TTTb re-enforces several facts that are important in 
attempting to understand where RMI phenomena take place. 
Instability vanishes in limit t p — > oo since the ratio of the radi- 
ation flux to energy density vanishes as well. The free energy 
for the instability is rooted in the anisotropy of the radiation 
field, which leads to a radiative flux that can then perform 
work on fluid fluctuations. Therefore, in the center or stars (or, 
e.g., at the surface of last scattering for the cosmic background 
radiation) RMI driving does not operate. Consequently, insta- 
bility favors regions relatively near the photosphere. 

The presence of v p h in the denominator of the left hand side 
implies that the instability favors magneto-acoustic oscilla- 
tions with small phase velocity. Given the oscillation wave- 
length, the rate at which radiation removes energy from a fluid 
disturbance is given by the phase velocity. That is, the rate at 
which radiative diffusion can remove energy from the wave is 
relatively small for relatively slow oscillations. 

3. NUMERICAL SETUP 

Our aim is to follow the nonlinear development of the 
RMI in a small patch of stellar interior. Below we describe 
the equations solved, numerical method, initial conditions, 
boundary conditions, and choice of parameters for simula- 
tions. 

3.1. Equations 

We solve the radiation-magnetohydrodynamic equations in 
the ideal MHD and grey diffusion limits 



dp 



(18) 



d 11 

+vVv= — V/7„+ (V xB)xB+— F + g(19) 

at p 4irp c 



3 kfP 

<9B 

— = Vx(vxB). 
at 



(22) 
(23) 



de 

— + v • Ve = ~7eV • \ + K a pc (E-aT*) 
at h 



(20) 



Our goal is to isolate the growth and nonlinear development 
of the RMI, which serves as justification for setting the flux- 
mean opacity Kp to a constant, as in the case of the Thomson 
scattering opacity. Radiation-hydrodynamic instabilities such 
as the ft -mechanism and the so-called strange mode instabil- 
ity (e.g.,(Glatzel 1994), which rely on opacity variations, are 
thus excluded. 

The Planck-mean and intensity-mean opacities normally 
enter independently into the radiation and gas energy equa- 
tions (eqns. lEOl - ED ). Here they are set equal to one another, 
implicitly assuming that the radiation f ollows a Planck dis- 
tribution with its own temperature (e.g., IBS03I) . We refer to 
these opacities collectively as the absorption opacity K a . In 
our calculations, the ratio K a / Kf is fixed to a value such that 
the condition for thermal locking (eq. 0) is satisfied. 

Equations dT8T>-(f23b are closed by an equation of state for 
the matter fluid, which we take to be an ideal gas of adi- 
abatic index 7. It follows that p g = pkT g /(pm p ) = (7- l)e, 
where p and m p are the mean molecular weight and proton 
mass, respectively. In the diffusion approximation, the radia- 
tion stress tensor P is related to the radiation energy density 
by V/E = 1/3, where 1 is the identity tensor. We denote the 
scalar radiation pressure by p T = E/3. No flux-limiters are 
employed (eq. 1221 ). 

A cartesian coordinate system is adopted, with the acceler- 
ation of gravity pointing in the negative z direction. Through- 
out this study, we restrict ourselves to two-dimensional simu- 
lations, with x denoting the coordinate direction transverse to 
gravity. 

3.2. Numerical Method 

To evolve the system of equations (fT8b-(f23]|, we em- 
ploy the publicly available finite-difference code Zeus-MP 
(lHayes et all 120061) . The MHD part is evolved using the 
Modified Charac t eristic s - Constrained Transport method of 
lHawlev & Stond d 19951) . Shock s are captured with the ar- 
tificial viscosity prescription of VonNeuma nn & Richtmveri 
( 1950). The radiation p art is evolved impl icitly via the Conju- 
gate Gradient method (Hayes et al. 2006), and interacts with 
the MHD sector via operator splitting. 

As a basic test o f our implementation, we verify that the 
initial conditions ( 33.31 ) do not evolve when unperturbed. 
We have also tested the elimination of the flux limiter in 
Zeus-MP by solving the diffusion equation in a unit square 
with periodic boundary conditions and constant initial E. A 
more stringent test involves comparing eigenfrequencies of 
magnetosonic eigenmodes on a uniform background in the 
therm ally-locked limit with the asymptotic values found by 
IBS03L Details on both of these tests are provided in Ap- 
pendixlBl 

3.3. Equilibrium Structure and Initial Conditions 

We initialize the computational domain with an atmosphere 
in hydrostatic, radiative, and thermal equilibrium. The equa- 
tions describing this atmosphere are obtained from equa- 
tions (TT9b-(f22b and the gas equation of state, by setting v = 0, 
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TABLE 2 

Models Evolved and Time- Averaged Properties 





Pm/Pg 


5R b 


e/e c b 




n x x n/ 


£kW(£r+£g) 


-Ek.v/Ek.h 


log 10 (Ap/p) f 


log 10 (AF/F) 


log 10 (AB/B) 


0.01 


0.3 


31 


1 .2 


2 90° 


128 


(2 ± 1.0)10 


i n-0 6 
10 


-3.0/-2.6/ - 1.7 


-3.0/-2.9/-2.0 


-2.0/ — 1.9/ -1.9 


0.03 






3.5 






(3±2.0)10~ 5 


10 -0.5 


-2.9/-2.5/-1.6 


-2.8/-2.9/-1.9 


—2.5/— 1.8/ — 1.8 


0.1 






11 






(7±3.0)10- 5 


10 -O.4 


-2.7/-2.3/-1.5 


-2.6/-2.7/-1.8 


-2.2/- 1.6/- 1.7 


0.3 






28 






(2±1.0)10" 4 


10 -0.2 


-2.1/-1.7/-1.3 


-2.2/-2.0/-1.5 


-1.5/- 1.3/- 1.5 


1 






60 






(3±2.0)10~ 3 


10 -0.4 


-1.3/-1.0/-0.8 


-1.7/- 1.4/- 1.0 


-0.7/-0.6/-0.9 


3 






91 






(1 ±0.6)10~ 2 


10 uo 


-0.7/-0.6/-0.5 


-1.0/- 1.0/-0.9 


+0.3/ -0.1/ -0.4 


10 8 






1 10 






> io- LV 










1 


0.01 


5.7 


9.8 


2 90° 


128 2 


(5±3.0)10~ 10 


10 u/ 


-4.6/-4.9/-4.2 


-4.4/-4.4/-4.4 


-4.5/-4.3/-3.5 




0.1 


18 


34 






(1 ±0.5)10~ 3 


10 -0.4 


-1.7/_ 1.4/- 1.2 


_1 9/_ 1.7/- 1.4 


-0.6/-0.6/-0.7 




1 


40 


79 






(5±2.0)10~ 3 


10 -0.8 


-1.0/-0.8/-0.6 


-1.3/-1.0/-0.9 


-0.9/- 1.0/- 1.2 




10 


40 


79 






(4±0.4)10 -3 


10 -3.0 


-0.8/-0.7/-0.7 


-1.4/- 1.3/- 1.5 


-1.8/-2.0/-2.4 


0.1 


0.3 


31 


11 


2 90° 


16 2 


(2±0.7)10 -6 


10 -0.6 


-3.0/-2.7/-2.1 


-3.1/-3.0/-2.5 


-Z6/-2.3/-2.3 












32 2 


(4±2.0)10" 5 


10 -0.4 


-2.7/-2.4/-1.5 


-2.8/-3.1/-1.9 


-2.3/- 1.8/- 1.8 












64 2 


(6±3.0)10" 5 


10 -0.4 


-2.7/-2.4/-1.5 


-2.7/-2.8/-1.8 


-2.2/- 1.7/- 1.7 


0.1 


0.3 


31 


11 


4 90° 


128 X 64 


(l.sio.sno- 4 


10 -0.6 


-2.3/-2.1/-1.4 


-2.6/-2.4/-1.7 


-1.9/- 1.5/- 1.5 










8 


256 x 64 


(i.o± 0.5)10-* 


10 -0.5 


-2.5/-2.2/-1.4 


-2.6/-2.6/-1.7 


-2.1/- 1.6/- 1.6 


1 


0.3 


31 


60 


2 75° 


128 2 


(7±3.0)10 -4 


10 -0.9 


-1.5/- 1.3/- 1.0 


-1.5/- 1.3/- 1.3 


—1.2/- 1.1/ — 1.2 










60° 




(2±1.0)10~ 3 


1Q -0.8 


-1.1/-1.0/-0.8 


—1.2/- 1.2/ — 1.1 


-0.9/-0.7/-0.8 



a Values of p Y j p*, p m / Pi, and e/e c at the center of the box. 

k The parameters and e c are determined by requiring that lo^ is 100 times the slow magnetosonic frequency of a mode with wavelength equal to // g , at domain center. 
c Width of the computational box in the direction transverse to gravity, in units of H g at z = 0. The vertical size is always 2//„. 

Angle between the imposed magnetic field and the background radiative flux. 
e Resolution in the directions transverse and parallel to gravity, respectively. 

^ Ratio of the root-mean-square fluctuation over the mean (eqns. |36| - |37| ) at positions zjL 2 = {0.2, 0.5, 0.8}, respectively. 
g The diffusion solver breaks down before saturation, due to order unity density fluctuations. 



: E, and dropping the time derivatives: 



dz 



dE 
"dl 
dlnp 
~dF 



= -3epg 

dlnpg 1 dln£ 
dz 4 dz 



(24) 
(25) 
(26) 



The equilibrium magnetic field B is taken to be uniform in 
space and therefore does not contribute to the equilibrium 
structure of the atmosphere. 

The problem is completely specified by nine parameters 
at one location. These are the ratio of box size L to the 
gas pressure scale height H g ; the Eddington ratio e; the ra- 
tio of radiation to gas pressure p r /p g ', the adiabatic index 7; 
the mean molecular weight /x, the ratio of isothermal sound 
speed to light speed c,/c; the optical depth over a gas pres- 
sure scale height t p ; the ratio of magnetic to gas pressure 
Pm/ Pg = B 2 / (&irp„); and the angle between magnetic field and 
radiative flux cos6» BF = B - F/(|B| |F|). 

We initialize the domain in such a way that there is cancel- 
lation to within a few parts in 10~ 6 between all the forces. This 
is possible because Zeus-MP computes gradients using finite 
differences on a staggered mesh. Evolving this initial state 
without additional perturbations excites vertical fast magne- 
tosonic modes, which decay (there is no RMI driving if the 
wave-vector is perpendicular to the magnetic field and/or par- 
allel to the radiative flux). 

Random initial perturbations in density are applied from 
z = Q.1L Z to z = 0.9L, for all x. The amplitude is 10~ 6 , and the 
random number is varied from cell to cell. The internal energy 



density is also perturbed following the density, Se/e = Sp/p, 
so that the gas temperature and radiation energy density re- 
main uniform everywhere. 

For simplicity, we consider atmospheres with £h = 1 , which 
are stably-stratified in the rapid diffusion limit (BS03). 

3.4. Boundary Conditions 

At the inner vertical boundary, we fix the density, gas-, and 
radiation energy density to their steady-state values, by con- 
tinuing the initial solution beyond the boundary. The radial 
velocity is reflected at this boundary, while the tangential ve- 
locity is set to zero in the ghost cells. The magnetic field is set 
to its uniform initial value. This arrangement ensures that all 
waves are reflected at the inner z boundary, and that the do- 
main does not collapse under the action of gravity. A constant 
radiative flux F is thus incident on the box from below. 

At the outer z boundary, we use an outflow boundary condi- 
tion that allows waves to leave the domain but avoids runaway 
mass loss. To account for the fact that this is an internal patch 
inside a hydrostatic structure rather than a domain with a 'vac- 
uum' external boundary, we set the outflow radial velocity so 
that the mass flux leaving the domain is that due do the ex- 
cess density over the steady-state value. At a given transverse 
coordinate, we thus set 



Vr.ghost : 



: max v 



w 0) 



max(p fl -po,0) 

Pa 



(27) 



where the subscripts and a label values in the initial and 
time-dependent solution at the last active cell inside the outer 
radial boundary. This prescription automatically includes the 
case where the density is much larger than the steady-state 
value, smoothly approaching a pure outflow boundary condi- 
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tion for the radial velocity. The density, gas energy density, 
magnetic field, and tangential velocity from the last active 
zone are copied into the two ghost zones. If the radiative flux 
between the last pair of active zones is positive, the radiation 
energy density in the ghost zones is set so as to keep this flux 
constant across cells. If instead the flux is zero or negative, 
then the radiation energy density is set to have zero gradient, 
yielding zero flux. 

The horizontal boundary condition is periodic for all vari- 
ables. 

3.5. Models Evolved 

We run two base sequences of models that vary the ra- 
tios of radiation to gas pressure p r / p g and magnetic to gas 
pressure p m /p g - These base sequences employ square com- 
putational boxes of two gas pressure scaleheights on a side 
(L x = L z = 2H p ), with a resolution of 128 x 128 cells, and 
a time step At equal to the diffusion time for a wavelength 
H g /S. The magnetic field points in the direction perpendicu- 
lar to gravity (i). This choice of field geometry allows waves 
to cross the domain many times while their amplitude grows. 
We perform a few runs that explore the effects of changing 
the box size, resolution, and field orientation. All models are 
shown in Table |2j with parameters showing values at the cen- 
ter of the domain. 

In all simulations, we take c/c,- = 1000, 7 = 5/3, and /j, = 1. 
The optical depth in the box is chosen so that the diffusion 
frequency c^diff is 100 times larger than the slow magnetosonic 
frequency for a wavelength equal to H g . The wave number de- 
pendence of cjditf ensures that all modes of shorter wavelength 
are well within the rapid-diffusion limit. This choice fixes the 
value of 5ft and thus e/e c . 

The ratio n- d / ftp is chosen so that 0% is 1000 times 
larger than the slow magnetosonic mode of wavelength equal 
to i/g/16. Our choice is motivated by results from tests 
of radiation-MHD waves on a uniform background (Ap- 
pendix |B). 

4. RESULTS 
4. 1 . Linear Growth and Saturation 

The general behavior of all models involves damping of the 
initially applied perturbations, and the subsequent emergence 
of unstable propagating waves. Because the background mag- 
netic field is perpendicular to gravity in most models, waves 
are allowed to cross the domain many timefl This leads to 
exponential growth of the kinetic energy. 

The evolution of the volume-integrated kinetic energy for 
a gas-dominated model {p r /p g = 0.1, p m / p g = 0.3, 128 2 ) is 
shown in Figure [3^, broken up into vertical and horizontal 
components, 



/ ( \pvt j clvd: 



E k .h= j ( -pv z x ) dxdz 



(28) 
(29) 



respectively. Initial density perturbations drive vertically- 
propagating waves that eventually decay. The horizontal ki- 
netic energy shows no such initial transients, growing out of 
numerical noise from the beginning of the simulation. Both 

4 To lowest order, the group velocity of slow magnetosonic modes lies 
along the magnetic field. 
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FIG. 3. — Panel (a): Evoluti on of the volume-integrated vertical and hori- 
zontal kinetic energies (eqns. 1281 - 1291 ) as a function of time, for the model 
with p r /p g = 0.1 and Pm/pg = 0.3 (0bf = 90°, resolution 128 2 ). Panel (b): 
rate of change of the horizontal kinetic energy in the same model, obtained 
by differencing variables at consecutive time steps (black). Also shown i s the 
sum of the terms that make up the kinetic ene r gy e quation (green; eq. 13 1 1 . 
Panel (c) shows each of these terms (eqns. 1321 - 1351 ') separately. 



energies saturate at nearly same time, closely following each 
other after the RMI overtakes the decaying waves. 
The equation governing the horizontal kinetic energy is 



d (\ 



dt\~2 pV > 



-V- 



1 2 



9 ( 

■^{Pg+Pr + Pn 



+-^(B-V)5 X . 



(30) 



By itegrating over the computational volume, we may write 



^k.h — ^adv ^sas ^sas Aor 3 



where 



4dv=" 
-^gas = " 
^rad = " 



v t ^dxdz 
ox 

v x —-dxd.z 

ox 

47r z \ dz dx 



(31) 

(32) 
(33) 
(34) 
(35) 



are the contributions from advection, gas pressure gradients, 
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radiation pressure gradients, and the Lorentz force, respec- 
tively. We combine the contribution from magnetic pressure 
and tension, since they cancel each other out to a large degree 
due to their phase offset during linear growth. 

Figure |3j) shows the evolution of the rate of change of the 
horizontal kinetic energy of the same model shown in panel 
(a). This rate is in excellent agreement with the sum of the 
terms comprising equation OTb . an additional verification of 
the accuracy of the solution method. 

The terms that make up the rate of change of the horizontal 
kinetic energy are shown separately in Figure [3}:. The gas 
and radiation pressure gradient terms grow at the same rate 
for ~ 15 e-folding times. The advective and Lorentz force 
components grow more slowly initially, but then increase their 
growth rate as the amplitude increases. 

Saturation in this model occurs when the Lorentz force 
reaches the same amplitude as the radiation pressure compo- 
nent. These two contributions are 90° out of phase, as ex- 
pected from the form of the RMI forcing in the linear regime 
(AppendixlAl; the saturation amplitude is small in this model 
(Sp/p ~ 3 x 10~ 2 ; Table ]2] >. This result is consistent with 
that of [Turner et all (120051) 7 who found that the RMI stopped 
growing when the field 'buckled'. The much lower radia- 
tive forcing in this model means that the amplitude of mag- 
netic field fluctuations required are much smaller. This bal- 
ance of stresses at saturation does not apply to other models 
with lower p T /p s , however. A more careful study of this in- 
terplay between magnetic and radiation stress tensors will be 
addressed in future work. 

The growth rates inferred from the evolution of are 
shown in Figure 0] for the sequence of models with vary- 
ing radiation pressure and fixed magnetic field {p m /Pg = 0.3; 
128 2 ; #bf = 90°). When radiation pressure is low, the mea- 
sured values agree within a factor of a few with the growth 
ra te of th e most unstable mode from the dispersion relation 
of BS03. Agreement improves as radiation makes an increas- 
ingly larger contribution to the total pressure, and the degree 
of forcing becomes larger. Comparing the growth rates mea- 
sured in the model with p T /Pg = 0.1 for resolutions of 64 2 and 
128 2 shows a decreasing agreement (at the few percent level) 
with coarser resolution, pointing to numerical dissipation as 
the likely cause of the discrepancy with the analytic value at 
low radiation pressure. We neverth eless take this as a confir- 
mation of the predictions of BS03 on the growth of the RMI 
when gas pressure dominates. 

Snapshots of the fractional change in density relative to the 
initial value during saturation are shown in Figure |5]for mod- 
els with fixed magnetic field and varying radiation to gas pres- 
sure ratios. When the radiation pressure is low, amplitudes 
are small and the flow is dominated by domain-filling struc- 
tures. The magnetic field becomes significantly deformed 
when p r > p m . A larger contribution from radiation pressure 
increases the saturation amplitude and amount of small-scale 
structure. 



The model with p, 



> 



10 and p m /pg = 0.3 deserves 



special attention. The amplitude of the fluctuations becomes 
large enough that the diffusion solver fails to converge. Future 
studies should address this region of parameter space using 
a suitable closure for the radiation moment equations and a 
larger computational domain. 

Another set of snapshots is shown in Figure [6] for models 
with fixed p r /p g and increasing magnetic field. Despite the 
fact that the contribution from radiation is significant, very 
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FIG. 4. — Linear growth rates inferred from the evolution of the horizontal 
kinetic energy as a function of the ratio of radiation to gas pressure, for the 
sequence of models with fixed magnetic field (squares). The solid line shows 
the growth rate of the most unstable mode from the dispersion relation of 
BS03. The dashed lines are the approxi mate growth rates of slow (red) and 
fast (blue) modes from equations (8) and jilt , respectively. 

small amplitudes are obtained when magnetic fields are low. 
In fact, no growth is seen at all for the lowest field config- 
uration when the resolution is 32 2 , in contrast to most other 
models. Decreasing the radiation pressure to p r /pg = 0.3 and 
keeping p m /Pg = 0.01 leads to no growth even at our base- 
line 128 2 resolution. Increasing the magnetic field strength 
increases the amplitude and spatial size of the dominant struc- 
tures up to the point where p m ~ p r . Further increases in mag- 
netic field cause the saturation amp litude to decrease sl ightly, 
giving rise to the structures seen bv lTurner etd](l2004b . 

4.2. Time-Averaged Properties 

We now discuss the quantitative behavior of the simulations 
in the saturated state. Saturation is defined as a period dur- 
ing which the kinetic energy is constant when averaged over 
timescales longer than the oscillation period or growth time. 
In some cases (e.g., p T /p g = 1), the system achieves an ini- 
tial phase of saturation lasting a few hundred sound crossing 
times, but then mass begins to leave the box at an increasingly 
rapid rate, changing its global properties. We always choose 
the time interval for averaging such that the mass in the box 
remains nearly constant. 

The time-average of the volume-integrated vertical kinetic- 
, horizontal kinetic-, gas-, and radiation energies is denoted 
by £k v , E^h, Eg, and E r , respectively. For field variables, we 
compute time- and horizontal averages 



1 



(TL X ) 



df / dxA, 



(36) 



where A(x,z,t) stands for a scalar such as density. The root- 
mean-square fluctuation is computed as 



AA: 



A 2 - (A) 



1/2 



(37) 



For all the models that reach saturation, we report in Table|2] 
the ratio of total time-average kinetic energyto time-averaged 
internal plus radiation energies E^ Aol /(E g + E r ), where Ek.tot = 
Ek,v + Ek,h, and the ratio of vertical to horizontal kinetic ener- 
gies Efy/E^f,. In addition, ratios of root-mean-square fluctu- 
ation to mean value of density Ap/p, vertical component of 
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FIG. 5. — Top: Instantaneous deviation of the density relative to its initial value during saturation, for models with increasing radiation pressure (as labeled) and 
constant magnetic field (pm/ Pg = 0.3, resolution 128 2 , #bf = 90°). Curves sho w a few magnetic field lines. Not e th at gravity points in the negative z direction. 
Bottom: Profiles of the time- and horizontal root-mean-square fluctuation (eq. 1371 ) divided by mean value (eq. 1361 ) of density (blue), radiative flux (red), and 
magnetic field (green), for the same models shown in the top panels. 




FIG. 6. — Same as Figure|5] but for a range of models spanning different ratios of magnetic pressure to gas pressure, as labeled. The radiation to gas pressure 
is fixed to p t / p g = 1, and the resolution is 128 2 . When the magnetic field dominates, the usual photon bubble behavior is recovered (Panel d). 
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FIG. 7. — Time-averaged properties as a function of the ratio of radiation 
to gas pressure, for models with Pm/Pg = 0.3 (#rf = 90°, resolution 128"). 
Panel (a) shows the ratio of total kinetic energy to the sum of the gas and 
radiation energies. Error bars show the size of the corresponding root mean 
square fluctuation. Panels (b), (c), and (d) show the ratio of the root mean 
square fluctuation to mean value at domain center of the density, vertical 
radiative flux, and magnetic field strength, respectively. 



the radiative flux AF/F, and magnetic field strength AB/B 
are provided at three points in the computational box. The 
lower panels in Figures[5]and[6]show that these averaged vari- 
ables change considerably with altitude, sometimes by up to 
an order of magnitude. To capture this variation, values at 
20%, 50%, and 80% of the box height are shown in Table 

Figure [7] shows the kinetic energy as a function p r /p g , for 
the sequence of runs with constant p m /Pg = 0.3 and reso- 
lution 128 2 . The kinetic energy is a monotonic function of 
the ratio of radiation to gas pressure, increasing the fastest at 
Pr/pg ~ 0.5. A similar trend is observed for the root-mean- 
square fluctuation of the density, radiative flux, and magnetic 
field strength, whose values at domain center are also shown 
in the same figure. Such increase in steady-state amplitudes 
is clearly the consequence of the fact that the radiation field 
is the ultimate source of energy ( §I2.31 >. The amplitude of the 
fluctuations in radiative flux and density are comparable to 
each other, while the magnetic field perturbations are larger 
by an order of magnitude at box center. This discrepancy is 
in part a consequence of the different spatial distributions of 
fluctuations in the box (e.g., lower panels of Figure|5}. 

The dependence of the saturated quantities on the ratio of 
magnetic to gas pressure is shown in Figure[8]for t he se quence 
of simulations with p t j p g = 1. As pointed out in §14.11 the ki- 
netic energy increases with field strength when the inequality 
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FIG. 8. — Time-averaged properties as a function of the ratio of magnetic to 
gas pressure, for models with Pi/p g = 1 (6*bf = 90°, resolution 128"). Panel 
(a) shows the ratio of the total kinetic energy to the sum of the gas and radia- 
tion energies. The ratio of the vertical to horizontal kinetic energies is shown 
in panel (b), and the normalized root mean square fluctuation at domain cen- 
ter of the magnetic field strength in panel (c). 



Pm -C Pi ~ p g obtains. Maximal amplitude is obtained when 
the contributions to the pressure are all comparable to each 
other. Further increase in the magnetic field suppresses ver- 
tical motions, leading to a lower total kinetic energy. This is 
clearly shown by the sharp decrease in the ratio of vertical to 
horizontal kinetic energy when going from p m /Pg= 1 to 10 
in FigurefSJ). The amplitude of the magnetic field fluctuations 
follow the magnitude of the kinetic energy. 

Table [2] also shows that varying the inclination of the mag- 
netic field relative to gravity changes the saturation amplitude. 
The optimal configuration appears to be a field perpendicular 
to gravity. When the magnetic field comes close to the ver- 
tical direction, the growth rates are eventually suppressed by 
the sine o f the angle between field and radiative flux (equa- 
tion |A20|. This change of saturat ion amplitude with angle 
was also seen bv lTurner et al.l (120041) . 

4.3. Dependence on Numerical Parameters 

We have selected a gas-dominated model (p r /Pg = 0.1, 
Pi/ Pg = 0-3, #bf = 90°) for studying the dependence of the 
saturation properties on the choice of numerical parameters. 
Figure [9] shows that the time-averaged kinetic energy in the 
box as a function of resolution is approaching convergence 
when using 64 cells per gas pressure scale height in each di- 
rection (128 2 cells in the box). The difference with using half 
the number of cells per direction decreases the kinetic energy 
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FIG. 9. — Time-averaged kinetic energy as a function of spatial resolution. 
Parameters of the model are p T / p g = 0.1, p m / Pg = 0.3, and 6>bf = 90°. Error 
bars denote the root mean square fluctuation. 

by ~ 15%. We thus consider our simulations are spatially well 
resolved. Note that using less than 10 cells per gas pressure 
scale height can significantly underestimate the saturation am- 
plitude; 

iTurner et al.1 ( 120041) found that the choice of time step is 
also relevant for correctly capturing the growth of the instabil- 
ity when using an implicit solver for the radiation sector. We 
have performed an additional set of simulations (not shown 
in Table [2]i that vary the time-step employed on the model 
with 64 2 resolution. No appreciable difference in the satura- 
tion amplitude is found when increasing the minimum time 
step from the fiducia l val ue, which resolves diffusion over a 
wavelength H g /16 ( §13 .5t . to the diffusion time for a wave- 
length H g /A . Minor differences of the order of 0.1 dex in the 
amplitude of fluctuations of field variables are the only differ- 
ence. One possible cause of this is that the nonlinear phase is 
dominated by large scale structures that are well resolved. 

The size of the box in the horizontal direction can change 
the value of the saturation amplitude by a factor of order unity. 
Table|2]shows that the simulation with 64 2 resolution has less 
than half the kinetic energy of that where the horizontal size 
is AH„ (128 x 64). This could be interpreted as a finite domain 
size effect, with some large scale waves being more capable 
of efficiently extract energy from the radiation field. Further 
enlarging of the box to 8// g in the horizontal direction causes 
a significant drop in the kinetic energy per unit volume, how- 
ever. We do not understand this effect at present. 

5. APPLICATION: MASSIVE STARS 

Massive main-sequence and blue supergiant stars have ra- 
diative envelopes, with an important contribution from radi- 
ation pressure to the overall support. Many OB stars have 
recently been found to have m agnetic fie lds up to kG strength 
or even higher (see, e.g.. iPetit et al.ll2012l for a recent compi- 
lation). Conditions near the surface of these stars are clearly 
susceptible to the RMlQ 

5 The analysis of BS03 assumes a constant magnetic field. The surface lay- 
ers of a star, where the instability is maximal, are such that the gas pressure 
scale height is much smaller than the radius. Stable magnetic fields in radia- 
tive envelopes are expect ed to appear mostly dipolar, w ith some contribution 
from higher multipoles (Braithwaite & Nordlund 2006). Such a field geom- 
etry varies on distances comparable to the stellar radius. Use of a constant 
magnetic field in local simulations is thus a good starting point to address the 



The winds of massive stars are known to be time-variable 
and show dumpiness over a wide range of spatial scales. For 
instance, mass loss rates inferred from UV resonance lines - 
which scale linearly with density - yield estimates that are at 
least an order of magnitude lower than those from collisional 
emission processes (e.g., Ha) , which depend on the density 
squared dFullerton et al.ll2006i) . 

Variability phenomena are traditionally divided into large- 
and small-scale. Below we elaborate on the potential impact 
of the RMI on each of these two categories. 

5.1. Large-Scale Variability 

It is currently thought that cyclic variability on scales com- 
parable to the stellar radius (e.g., time-dependent absorption 
superimposed on the blue edge of P Cygni profiles) origi- 
nates in so-called Corotating Interacting Re gions (CIRs), in 
analogy with the Solar Wind (Mu llanlll984l) . In this model, 
azimuthal fluctuations at the base of the wind cause faster 
streams to overtake slower ones, generating shocks and vari- 
ations in the optical depth that follow a spiral pattern due 
to angular-momentu m conservation (e.g. JCranmer & Owockil 
1996, Dessardl200l . 

The source of these fluctuations in massive stars is not 
established yet. Proposed mechanisms include flow in a 
spatia lly inhomogeneous magnetic field (lUnderhill & Fahevl 
1984) and multi-periodic non-radial pulsations ( Kauf eret al.l 
2006), the latter possibly e xcited by near surfac e convection 
zones from opacity peaks (Cantiel lo et al.l 12009 ). In at least 
one source (HD 64760), the spectral line variability can be 
accou nted for if the spots do n ot corotate with the stellar sur- 
face (lLobel & Blonimel[2()08 ). lending support to the multi- 
periodic non-radial pulsation hypothesis. From our results, 
the RMI is a clear candidate as a source of such pulsations. 

5.2. Small-Scale Variability 

Line-driven winds are subject to the L ine Deshadowing 
Instab i lity (LDI) above t he ph o tosphere (|Lucv & Solomonl 
fl970l iMacGregor etaf] [T979t lOwocki & Rybickil 1 19841) 7 
When the velocity of a parcel of gas in the wind is perturbed 
outward, the corresponding Doppler shift moves the relevant 
rest-frame transitions into an energy range where the contin- 
uum flux from the star is higher, increasing the outward accel- 
eration. The instability tends to favor small spatial scales. The 
nonlinear development of the LDI results in extended regions 
of sharp density and velocity fluctuations that le ad to the com- 
pressi on of wind material into thin radial shells (lOwocki et al.l 
1988). These shells are subsequently brok en up by secondary 
instabilities, resulting in a clumped wind (iDessart & Owockil 
2003). It is currently thought that the LDI can account for 
most of the phenomenology requiring small scale inhomo- 
geneities, such as transient emission line substructures, black 
troughs of saturated UV resonance l ines, discrepan t mass-loss 
estimates, and X-ray emission (e.g.. lFuUertonll20TTl) . 

A remaining question is the presence of clumping at the 
base of the wind. Simulations of the 'self -excited' LD I find 
that damping due to scattered radiation dLucyl 1984) sup- 
presses the onset of fluct uations below ~ 1.5 stellar radii 
dRunacres & Ow ocki 2002). A number of observational stud- 
ies have found evidence f or clumping much cl oser to the base 
of the wind. In particular, Nai arro et al.l (1201 ll) used different 
line diagnostics to reconstruct the radial profile of the clump- 
ing factor (average of p 2 over the squared average of p) for 

effects of the RMI in stars with radiative envelopes. 
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the entire wind inside ^100 stellar radii in £ Pup, finding 
that th e peak of this factor occurs inside 2 stellar radii. Re- 
cently, Sundqv ist & Owoc ki (2012) have repo rted LDI simu- 
lations that agree with the Naj arro et alJ (1201 ll) clumping pro- 
file when including limb-darkening effects (which decrease 
line-drag from scattered radiation) and generic sound wave 
perturbations at the base of the wind. 

The RMI can readily drive sub-photospheric density pertur- 
bations on small spatial scales. IBS03I found that in the ther- 
mally locked limit (eq. (3)) the instability operates up to a 
wave number 



(CVhS/c? 
(C 5 ^thg/v 



1/2 



(slow) 
(fast) 



(38) 



where £ is given by equation (0. For a near equipartition 
magnetic field, this translates into a minimum unstable wave- 
length of 



A cut ^6(p_^3)- 1/2 ^ 6 2 (^ 



1/2 



1/2 



km, (39) 



K es 

where p_g is the density in units of 10~ 8 g cm" 3 , #3 the accel- 
eration of gravity in units of 10 3 cm s" 2 , c, 6 the isothermal 
sound speed in units of 10 km s" 1 , and K es is the electron scat- 
tering opacity. Comparing to the gas pressure scale height at 
the photosphere, 

-1 



#g,ph ~ («FPph) 

~3000pl8 (^) km, 



(40) 
(41) 



shows that the RMI in its WKB version operates over a dy- 
namic range of at least 100 in wavelength. 

6. SUMMARY 

In this paper we have studied the nonlinear development of 
the Radiation-Driven Magneto- Acoustic Instability (RMI) in 
the regime where the gas pressure is comparable or dominates 
over that due to the radiation and magnetic field, and in the 
limit of rapid heat exchange between matter and radiation 
(eq. J3|). Our primary results are: 

1 . - We have confirmed the predictions of Blaes and Socrates 
(2003). In particular, the RMI instability operates when gas 
pressure dominates over radiation pressure and when the 
magnetic field is weak. The growth rate of the kinetic energy 
is consistent with that expected for the most unstable mode 
(Fig.©. 

2. - The saturation amplitude is an increasing function of 
the ratio of radiation to gas pressure (Fig. |7]l over all of the 
dynamic range studied here, p T / p g £ [10~ 2 , 10]. At the lowest 
end of this interval, the saturation amplitude is small but 
non-negligible (5p/ p ~ 10~ 3 ). 

3. - The saturated state varies non-monotonically with 
magnetic field strength (Fig. [8}. For moderate forcing by 
radiation (p r ~ p„), the saturation amplitude increases with 
field strength as long as magnetic pressure is sub-dominant. 
Once the magnetic field dominates, vertical motions are 
suppressed and increasingly smaller amplitudes are obtained. 



4. - The RMI may be a source of sub-photospheric per- 
turbations in the envelopes of massive main-sequence and 
blue supergiant stars. The instability can serve as a source 
of large spatial scale perturbations for Corotating Interacting 
Regions. The RMI also operates over a wide dynamic range 
in wavelength, providing sufficient seed perturbations for the 
Line Deshadowing Instability even if a convection zone is 
absent. 

The computational requirements to capture the RMI in real- 
istic contexts are very demanding. The instability is strongest 
in regions near the photosphere, where diffusion is fast rela- 
tive to the mode period and where the pressure scale height is 
much smaller than the stellar radius. Capturing the instabil- 
ity requires resolving diffusion over the pressure scale height 
in space and in time, introducing an additional computational 
cost. For example, the ratio of the sound crossing time over 
the diffusion time at a given wavelength is pro port ional to the 
diffusion speed over the sound speed (oc 5ft, eq. IT4l ). For mas- 
sive stars this ratio can be as high as ~ 10 5 , requiring a com- 
parable number of diffusion steps per Courant time to capture 
the RMI. 

Future work will address the saturation of the RMI, and fur- 
ther application to astrophysical systems. 
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APPENDIX 

A. DERIVATION OF LINEAR GROWTH RATES AND 
STABILITY CRITERIA 

A. 1 . Basic Magneto-Acoustic Waves 

In the absence of stratification and a changing radiation 
field, compressive fluid perturbations are forced by a com- 
bination of their own acoustic response and the Lorentz force. 
In this limit, conservation of momentum to linear order can 
be written as 



-iujQp5\ = — /k 



5p g + - 



B <5B 



An 



+ i™5B 

4n 



(Al) 



where 6 denotes Eulerian perturbation, and Sp g = c 2 5p when 
ST — > (eq. 0). We assume that all perturbed quantities 
take plane-wave form and are oc exp (k • x - ojt ). Using the lin- 
earized mass conservation equation, the linearized induction 
equation 

B . „ . (k B) 



<5B = — (k-Sv)- 



CJ 



■6y, 



(A2) 



and the component of the linear ized momentum equation 
along B, one can rewrite equation (lAlb as 



UJQ p 



Cj 2 



—^cf + vi k-(k-v A ) v A 



(A3) 
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where ui 1 = oji - (k ■ v A ) 2 . Projecting onto k and using mass 
conservation yields the magnetosonic dispersion relation, 



*(K-F) — -Wdifll — 



5T, (A15) 



2 — w 



"ph- 



1 

2 



(c 2 + vi)± 



(cf + 



/ 2 ) 2 -4(/c-v A ) 2 c 2 



(A4) 
(A5) 



A. 2. Radiative Correction: Driving and Diffusive Damping 

We now include stratification in the background state and 
allow the radiation field (or equivalently the temperature) to 
vary. The exchange of heat between radiation and energy is 
treated as a correction to the mode frequency. We thus have 



uj = wq + zT, 



(A6) 



where T is a driving or damping rate such that T <C uiq. The 
linearized momentum equation is now 



-ik 



5p g + 6p r + 



B SB 

4n 



k B 

+ i——SB+g5p, 

4-7T 



(A7) 



where £ = i8v /uj is the Lagrangian displacement, and p x = 
E/3 = aT 4 /3 is the radiation pressure in the optically thick 
limit. 

Defining P = p„ + p T , using the condition of hydrostatic 
equilibrium pg = VP, and expressing Eulerian perturbations 
in terms of the Lagrangian displacement 



5T 



<5B = /(k-B)|-/(k-|)B, 
one can re-express momentum balance as 



-uj 2 + cj ( 2 ) | = 5T+ i c 2 £ x (k x V In p) 



where 



dT 



tk5r+(k-ovr] 



(A8) 
(A9) 
(A10) 

(All) 
(A12) 



contains the driving and dam ping t erms. The second term on 
the left hand side of equation (I A 1 11 > includes the terms arising 
from acoustic and magnetic forces (equation lA3l l. The second 
term on the right hand side is perpendicular to the Lagrangian 
displacement and hence cannot do work on the perturbation 
(e.g., [Socrates et al. 2005). Terms of higher order in {kH g )~ l 
have been discarded. 
To first order in V / u>o, the growth rate is 

2 6\ ■ S\ (3) 

Before evaluating equation dA13b . it is instructive to ex plicitly 
separate the driving and damping terms in equation (1A121 >. 
To do so we solve for the temperature perturbation in the 
linearized equation for the total energy density. Defining 
U = e + E, we have 



- iui6U + S\ ■ VU+i(U+P)k ■ S\ = -ik ■ 5F 



(A14) 



or 



+ nTS\ ■ Vs = 



where 
T (ds 



dU 
dp 



U+P 1 / 4 

p p V 3 



(A16) 



To first order in uj/ujam and in (kH g ) the temperature per- 
turbation is 



5T_ 
T 



uj nT 


( ds \ 


+ k.Fl 


Wdiff E 


\d\npj 


j ^diffP . 



bp 

p 



— ■ (A17) 



The first term represents damping by radiative diffusion 
(eq. 1A16I ): entropy is lost upon compressioifl Note that 
we h ave d ropped the first term on the left hand side of equa- 
tion ( IA15b : it becomes important if (p T /pg) < (w/wdiff), con- 
tributing with a global factor of order unity to the growth rate 
(e.g., compare r s i ow and the maximum growth rate from the 
dispersion relation in Figure HI the choice of parameters is 

such that wo/wdiff = 0.01). 

We now evaluate equ ation (IA13b by inserting equa- 
tion ( IA17I) in to equation ( IA12I ). projecting onto 5v, using 
equation iA3i to eliminate a term proportional to (k x S\), 
and keeping the same order of approximation, we find 



Sp\ 2 (8P/dT) p f uj 
p J (dE/dT) p \uim 



ui x 



nT 



ds 
dlnp 



(k-v A )-(kxv A )-(kxF) 



(A18) 



The second term inside the square bracke ts represents driving 
by the background radiative flux (BS03). Note that only the 
component of the flux perpendicular to k contributes to the 
driving, and that driving vanishes when the magnetic field is 
completely par allel or perpendicular to the wave v ector. By 
using equation ( IA3t . the denominator of equation dA13l > can 
be written as 



S\ ■ Sv = 



•vi) 



(A19) 



Equations ©, ( lA13l . dATBt . (lA16l i and (KT9[ then yield the 
growth rate to first order in (T/u>o), (wo/wdiff), and {kH g )~ l , 



T~ - 
2 



" P h 



■(t-v A ) 



2v 2 - 



3 Pg 
4E 



+ 1 



CVph 



(k-v A )(kxv A )(kxF) 



[v 2 h - 



(fc-v A ) 2 ] 



Vph 



(A20) 



Equation ( IA20I ) is identical to the result obtained by lBS03l in 
the limit of rapid heat exchange (their equation 93). 

The stability criterion is o btaine d by requiring that the last 
squared bracket in equation ( IA20b is positive. Ignoring angu- 
lar factors, we can write 



F>5[w,k,B] h> g +-£l vm, (A21) 

6 In cosmology, this effect is known as "Silk damping", after Silk 1 1968). 
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FIG. 10. — Deviation from unity of the normal ized projection a n of the 
numerical solution to the diffusion equation (eq. IB31 ). evaluated at ; = to- 
The left panel shows results for different spatial resolutions, using a time step 
equal to the diffusion time at the grid scale. The right panel shows results 
with 128 cells per side of the square domain, using a minimum time step 
equal or larger than the diffusion time step. 



where 



S[w,k,B] 



-(k-y A ) 2 



(A22) 



measures the ratio of compressional energy to kinetic energy 
due to motion along B. 

B. TESTS OFZEUS-MP 

The magnetohydrodyn amic capa bilities of Zeus-MP have 
been extensively tested (Hayes et al. 2006). However, the ra- 
diation module has only been tested when coupled to the hy- 
drodynamics solver. Here we describe additional tests that 
quantify the accuracy with which equations (TT~8T>-(f23l> are 
solved. 

B.l. Diffusion Solver 

To test the removal of the default flux limiter in Zeus-MP, 
we compare the analytic and numeric evolution of a sinusoidal 
thermal wave. In this test all the accelerations, magnetic 
fields, and absorption opacities are suppressed in the code. If 
the density and gas pressure are set initially to uniform values, 
and the velocities vanish, the only time-dependent equation 
left is the diffusion of radiation energy density: 



dE 

— -V-(DV£) = 0, 
at 



(Bl) 



where D = c /(3kpp) is the dif f usion coefficient. 

Following iTurner & Stond (|2001), we employ a periodic 
square box of length L and a constant diffusion coefficient. 
A sinusoidal thermal wave with wavelength equal to the side 
of the box should evolve according to 

E a (x,y,t) = Eo + Ei e-'l'° sin (lirj) sin (2tt|) , (B2) 

where Eq and E\ are free parameters, the subscript a denotes 
analytic solution, and to = L 2 /(&ir 2 D) is the e-folding time of 
the solution. At t = 0, we in itialize the radiation energy density 
according to equation (IB21 >. 



To diagnose the agreement between the numeric and an- 
alytic solution, the numerical radiation energy density E n at 
t = to is projected onto the appropriate Fourier mode and nor- 
malized 



4 e 1 
Z2£7 



JO 



E n (x,y,to) sin (Zttj^J sin f^TT—J dxd;y. 



(B3) 

The quantity a n - 1 is a measure of the global deviation of the 
numerical solution from the analytic expectation. 

A number of tests were performed using different spatial 
resolutions and time steps, with an initial amplitude E\ /Eq = 
0.5. Results are shown in Figure [TO] Agreement is better than 
1% for the chosen set of parameters. 

B.2. Thermally '-Locked Radiation-MHD Waves on a 
Uniform Background 

Here we test whether the solution of equations (fT~8b-(f23T> 
by Zeus-MP correctly attains the limit of rapid diffusion and 
rapid heat exchange (eqns. 13 and (3)) in the presence of mag- 
netic field fluctuations. To this end, we follow the evolution 
of an individual fast magnetosonic eigenmode in a uniform 
square d omain of size L w ith periodic boundary conditions 
(see also iDavis et ail 120 121 and Jiang e t al J 120121 for similar 
tests). 

In the absence of stratification, radi ation diffusion causes 
damping, with an asymptotic rate (eq. MA20I ) 



r ■ =- 

x urn 



1 



(!) 



2v 2 - 

ph 



'ph -(fc-VA) (Pg_ + l \ ( 1 + fM £» 
4 Pr J\ pj V 
(B4) 

where Tb ox = k?pqL is the optical depth in the box. The rapid 
diffusion and heat exchange conditions ensure that the real 
part of the frequency ujq is that of an is othe rmal magnetosonic 
mode in the absence of radiation (eq. 0A41 ). 

The parameters of the problem are a wavelength equal to 
the domain size, wave-vector along the z axis, and a magnetic 
field at 45 degrees to this axis. The fiducial model has a ra- 
tio of pressures p T /pg = 10~ 2 , p m /pg = 0.1, an optical depth 
such that Wdifr/wo = 100, a ratio of absorption to flux-mean 
opacity such that Wth/wrj = 10 5 , a resolution of 32 cells per 
wavelength in both directions, and a time step equal to the 
diffusion time at the grid scale. All other parameters have the 
default values described in 33.51 The domain is initialized 
with a plane wave in density and gas pressure with amplitude 
Sp/po = Sp/po = 10~ 2 , and spatial dependence oc cos(27rz/L). 
The initial magnetic fi eld and v eloc ity perturbations are ob- 
tained from equations ( IA2b and ( lA31 l. respectively. 

To quantify agreement with the analytic expectation, we 
project the fractional deviation of the density from its equilib- 
rium value onto the appropriate Fourier component, obtaining 
a time-dependent coefficient 



2 

u 



o Jo 



l r 



p (x,z,t) 

P0 



1 



cos( ^ ) clvd;:. (B5) 



The above coefficient is then fitted with a cosine function of 
exponentially decreasing amplitude, yielding an oscillation 
frequency, damping rate, and initial amplitude. 

We evolve two series of runs. The first has no magnetic 
field, and samples different spatial resolutions, time steps, and 
absorption opacities. Results are shown in Figure[TTh-c. in the 
form of a fractional difference between the analytic and fitted 
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FIG. 11. — Properties of rapidly-diffusing, thermally-locked radiation-magnetosonic modes on a uniform background, as a function of system parameters. 
Panels (a)-(c) show the fractional difference between the fitted and analytic eigenfrequencies (upper and lower panels for damping rate T and oscillation frequency 
wo, respectively) obtained when varying the degree of thermal locking and numeric parameters. The baseline model has p r /pg = 10 , no magnetic field (va = 0), 
ui L \iff/uiQ = 100, uifa/ujQ = 10 5 , A/Ajc = 32, and timestep equal to the diffusion time at the grid (At = Afjjff). Panels (d)-(e) show eigenfrequencies as a function of 
the relative strength of gas, radiation, and magn etic pre ssure . The baseline model is the same as for panels (a)-(c) but with a non-zero magnetic field (pm/ p% = 0.1 
in panel d). Curves show analytic values (eqns. lB4l and lA4t while symbols show results from fits to simulations. 



eigenfrequencies. Agreement is better than 1 % in ojq and 10% 
in r un i, respectively, for all chosen parameters. The second set 
has a non-zero magnetic field, and samples different radiation, 
magnetic, and gas pressures. Results are shown in FigurefTTtl- 



e. Agreement between analytic and fitted eigenfrequencies is 
excellent over a wide region of parameter space. The only 
exception is the damping rate for very large magnetic field, 
where a larger damping rate is measured. 



REFERENCES 



Arons, J. 1992, ApJ, 388, 561 

Begelman, M. C. 2001, ApJ, 551, 897 

Blaes, O., & Socrates, A. 2001, ApJ, 553, 987 

Blaes, O., & Socrates, A. 2003, ApJ, 596, 509 

Braithwaite, J., & Nordlund, A. 2006, A&A, 450, 1077 

Cantiello, M., Langer, N., Brott, I., de Koter, A., Shore, S. N., Vink, J. S., 

Voegler, A., Lennon, D. J., & Yoon, S.-C. 2009, A&A, 499, 279 
Catlett, C. e. a. 2007, in HPC and Grids in Action, ed. L. Grandinetti 

(Amsterdam: IOS Press) 
Cranmer, S. R., & Owocki, S. P. 1996, ApJ, 462, 469 
Davis, S., Blaes, O., Turner, N., & Socrates, A. 2004, in ASP Conference 

Series, Vol. 311, AGN Physics with the Sloan Digital Sky Survey, ed. 

G. T. Richards & P. B. Hall, 135 
Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJS, 199, 9 
Dessart, L. 2004, A&A, 423, 693 
Dessart, L., & Owocki, S. P. 2003, A&A, 406, LI 

Fullerton, A. W. 2011, in IAU Symposium, Vol. 272, IAU Symposium, ed. 

C. Neiner, G. Wade, G. Meynet, & G. Peters, 136-147 
Fullerton, A. W., Massa, D. L., & Prinja, R. K. 2006, ApJ, 637, 1025 
Gammie, C. F. 1998, MNRAS, 297, 929 
Glatzel, W. 1994, MNRAS, 271, 66 

Hawley, J. F, & Stone, J. M. 1995, Comp. Phys. Comm., 89, 127 

Hayes, J. C, Norman, M. L., Fiedler, R. A., Bordner, J. O., Li, P. S., Clark, 

S. E., ud Doula, A., & Low, M.-M. M. 2006, ApJS, 165, 188 
Hsu, J. J. L., Arons, J., & Klein, R. I. 1997, ApJ, 478, 663 
Jiang, Y.-F, Stone, J. M., & Davis, S. W. 2012, ApJS, 199, 14 
Kaufer, A., Stahl, O., Prinja, R. K., & Witherick, D. 2006, A&A, 447, 325 
Klein, R. I., Arons, J., Jernigan, G., & Hsu, J. J.-L. 1996, ApJ, 457, L85 
Lobel, A., & Blomme, R. 2008, ApJ, 678, 408 



Lucy, L. B. 1984, ApJ, 284, 351 

Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879 

MacGregor, K. B., Hartmann, L., & Raymond, J. C. 1979, ApJ, 231, 514 

Mullan, D. J. 1984, ApJ, 283, 303 

Najarro, F, Hanson, M. M., & Puis, J. 2011, A&A, 535, A32 
Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 
Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 

Petit, V., Owocki, S. P., Wade, G. A., Cohen, D. H, Sundqvist, J. O., Gagne, 
M., Mai'z Apellaniz, J., Oksala, M. E., Bohlender, D. A., Rivinius, T., 
Henrichs, H. F, Alecian, E., Townsend, R. H. D., ud-Doula, A., & the 
MiMeS Collaboration. 2012, MNRAS, in press, arXiv:121 1.0282 

Prendergast, K. H, & Spiegel, E. A. 1973, Comments on Astrophysics and 
Space Physics, 5, 43 

Runacres, M. C, & Owocki, S. P. 2002, A&A, 381, 1015 

Silk, J. 1968, ApJ, 151,459 

Socrates, A., Blaes, O., Hungerford, A., & Fryer, C. L. 2005, ApJ, 632, 531 

Socrates, A., Parrish, I. J., & Stone, J. M. 2008, ApJ, 675, 357 

Sundqvist, J. O., & Owocki, S. P. 2012, MNRAS, in press, arXiv:1210.1861 

Tao, T., & Blaes, O. 201 1, ApJ, 742, 8 

Turner, N. J., Blaes, O. M., Socrates, A., Begelman, M. C, & Davis, S. W. 

2005, ApJ, 624, 267 
Turner, N. J., Quataert, E., & Yorke, H. W. 2007, ApJ, 662, 1052 
Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 

Turner, N. J., Yorke, H. W., Socrates, A., & Blaes, O. M. 2004, RvMxAA, 
22, 54 

Underhill, A. B., & Fahey, R. P. 1984, ApJ, 280, 712 

VonNeumann, J., & Richtmyer, R. D. 1950, Journal of Applied Physics, 21, 
232 

Wade, G. A., et al. 2012, MNRAS, 425, 1278 



